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\Q Abstract 

We study the internal resonance, energy transfer, activation mecha- 

1 | nism, and control of a model of DNA division via parametric resonance. 

While the system is robust to noise, this study shows that it is sensitive 
to specific fine scale modes and frequencies that could be targeted by low 
intensity electro-magnetic fields for triggering and controlling the division, 
i-^ The DNA model is a chain of pendula in a Morse potential. While the 

(possibly parametrically excited) system has a large number of degrees 
of freedom and a large number of intrinsic time scales, global and slow 
i i variables can be identified by (i) first reducing its dynamic to two modes 

exchanging energy between each other and (ii) averaging the dynamic of 
t-H the reduced system with respect to the phase of the fastest mode. Sur- 

prisingly the global and slow dynamic of the system remains Hamiltonian 
\l (despite the parametric excitation) and the study of its associated effec- 

tive potential shows how parametric excitation can turn the unstable open 
state into a stable one. Numerical experiments support the accuracy of 
the time-averaged reduced Hamiltonian in capturing the global and slow 
dynamic of the full system. 

{vq In this paper we study the internal resonance, energy transfer, activation 

y—i mechanism, and control of a model of DNA division via parametric resonance. 

While DNA macro-molecules are robust to noise, our study shows that they 
• *-j are sensitive to specific fine scale modes and frequencies that could be targeted 

K% by low intensity electro-magnetic fields for triggering and controlling the divi- 

sion. The suggested method of control is supported not only by the observation 
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that DNA vibrations induced by electric-fields or microwave absorbtion are an 
experimental reality but also by the fact that electric field-induced molecular 
vibrations have already been used as a noninvasive cell transfection protocol. 
Our study also raises the question on whether enzymes are using the proposed 
mechanism to initiate the opening of DNA strands. This question is to put 
into correspondence with increasing theoretical and experimental evidence that 
low-frequency vibrations do exist and play significant biological functions in 
proteins, DNA molecules, and other bio-macromolecules. 

1 Introduction 

The model considered in this paper is a chain of pendula in a Morse potential, 
with torsional springs between pendula, that mimic real DNA characteristics 
[24| [38] . Previous studies EH El |9] , mainly numerical, showed that this model 
exhibits an intriguing phenomenon of structured activations observed in many 
bio-molecules: while the system is robust to noise, it is sensitive to certain 
specific fine scale modes that can trigger the division. Below, we will describe 
briefly the results of our analytical study on this intriguing phenomenon and 
our effort in the control of this DNA model via parametric resonance. 

By using the Fourier modal coordinates [H E2] , this model can be seen as 
a small nonlinear perturbation of n harmonic oscillators with frequencies uj a — 
-y/2(l — cos 2-Ka/n), a — Q,...,n— 1. The coarse variable of the (approximate) 
Oth mode, which is the average angle of the pendula, corresponds to the angle of 
the frame defined by the radius of gyration in our molecular studies (40J ISJ S2] ■ 
This variable is a global and slow variable and it plays an important role as the 
reactive coordinate for our DNA work. Moreover, this reactive mode forms a 
nearly : 1 resonance with any other mode, each of which has an O(l) frequency. 
This fact leads to small denominators or coupling terms in the corresponding 
averaged equations or normal forms [26l [27l ES ESI [13 H31 EH |20l [2]. Since 
other modal frequencies are not rationally commensurate or have significant 
time scale separation, we do not expect strong resonance among them. Extensive 
simulations confirmed our expectation. We observed: (i) the energy transfers 
mainly from an excited mode to the reactive mode, triggering the division, (ii) 
only an extremely small amount of energy transfers from the excited mode to 
one or two other modes via near resonance. This observation, together with 
a rigorous error estimate jS], show that two-mode or three-mode truncation, 
i.e., a truncated system that includes the reactive mode, the excited mode, and 
perhaps a mildly affected mode should provide an adequate reduced model for 
our analytical study on the activation mechanism. 

By applying the method of partial averaging [El [fl ESI ESI ESI 1131 E5] for 
nearly : 1 resonance, we obtained the averaged equations for the reduced mod- 
els of this chain of Morse oscillators up to nonlinear terms of very high degree. 
These averaged reduced equations not only reveal the coupling between the ac- 
tion of the excited modes and the dynamics of the reactive mode, they also shed 
lights on the phase space structure of the activation mechanism. Moreover, they 
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allow us to estimate analytically the minimum activation energy for each excited 
mode. These analytical estimates not only match well with those obtained from 
simulations of the full DNA model, but also uncover a relationship between 
the frequency of the excited mode and its corresponding minimum activation 
energy. These findings also provide an analytical and deeper understanding of 
the internal mechanism that is responsible for the phenomenon of structured 
activations. 

Based on our understanding of its internal dynamics, we introduce a method 
for controlling the division of this DNA model via parametric excitation, that is 
in resonance with its internal trigger modes. By choosing appropriate external 
excitations and frictions, we uncover a class of trajectories that show how the 
parametric resonance can be used to drive the averaged reduced model from 
its (almost) equilibrium state to its open state. The identification of an effec- 
tive Hamiltonian (after reduction and partial averaging) opens the possibility of 
studying the global phase space structures of our averaged reduced model and 
sheds lights on the significant trajectories mentioned above. Moreover, we ex- 
tend the results for the averaged reduced model with parametric resonance and 
friction to the reduced as well as the full model with parametric resonance and 
friction. These findings support the conjectures that (i) low intensity electro- 
magnetic fields can be used with parametric resonance to inject energy into the 
trigger modes and destabilize the DNA chain and keep it near the open state 
for replication and transcription, (ii) enzymes may use similar method to initi- 
ate open state dynamics for DNA replication and transcription. Although the 
issues of inhomogeneity, helicity, and environmental effects (such as noises) will 
be ignored for now, they will be taken into consideration in our future work. 

2 A Model of DNA Division and Structured Ac- 
tivations 

Our analytical study has been inspired by the work of Mezic and Eisenhower 
at UC Santa Barbara and Marsden and Du Toit at Caltech. Their studies [231 
8, 9,[TT], mainly numerical, showed that this DNA model exhibits an intriguing 
phenomenon of structured activations. While Eisenhower and Mezic did apply 
Arnold's method of partial averaging [l] to a truncated Hamiltonian in their 
analytical study of a chain of Duffing oscillators [101 HU H2] , they did not extend 
it to the DNA model. The reason that they gave in [T3j for studying the Duffing 
case was that "The exponential form of the Morse potential makes analytical 
progress difficult,..." Since we have not seen this kind of simplification in the 
DNA literature, we will, in this study keep the Morse potential (and the added 
difficulty) . 

2.1 A Model of DNA Division 

The model was first introduced in [M] . It is a chain of equivalent pendula that 
rotate about the axis of a fixed backbone with an angle Ok measured from the 
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Fi gure 1: (a) A model of DNA division with 30 pendula. (b) The phase space (0,po) of 
a single pendulum in a Morse potential without coupling. It has two stable equilibria and a 
saddle. The black curve is the homoclinic orbit that separates two types of motion, namely, 
the libration near the equilibria (# e ,0) = (±2.346,0), and the flipping across the saddle (0,0). 



upward position. The pendula interact with nearest neighbors along the back- 
bone through harmonic torsional coupling, and with pendula on the opposing 
immobilized strand through a Morse potential that has two stable equilibria and 
a saddle. See Figure [I] 

The Hamiltonian that describes the motion of these n coupled pendula is 
given by 



k=l 



z% + ls{6 k - fc _O 2 + ( e -«-Wi+co.fc)-*o} _ i 
2mn z 2 2adh \ 



(1) 

with periodic boundary condition, 6 = 9 n . The first term is the kinetic en- 
ergy terms of n-pendula. The second term is the torsional coupling terms. The 
third term is the Morse potential terms, which model the hydrogen bonds of 
the respective DNA base pairs. In Eq. ([I]), m and h represent the mass and 
the length of each pendulum, and Pgk = mh 2 (d9k/dt) is the generalized mo- 
menta conjugate to df.- The parameter S determines the strength of the nearest 
neighbor coupling, while the parameter D determines the strength of the Morse 
potential. The parameter xq determines the equilibrium distance of the Morse 
potential, while the parameter ad determines the width of the Morse potential. 
All the parameter values were chosen to best represent the typical values for the 
opening and closing dynamics of DNA division [H [22] as follows: m = 300 amu 
(typical mass of a DNA base (nucleotide)), h — 1 nm (typical radius of DNA), 
S = 42 eV, D = 0.42 eV, x = 0.3 nm, and a d = 7 nm" 1 . 

After dividing the both sides of Eq. ([!]) by S, the Hamiltonian can be non- 
dimensionalizcd as 



fc=i 



mPe) = E + 5<°k- ° k -if + e (e-a-H-'-.-*) - l)' 



(2) 



where pek = dO^/dr is the dimcnsionless momentum defined with respect to 
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the dimensionless time r = J S/mh 2 t. Thus, in the present study, one unit 
time (t = 1) corresponds to t = yfmh 2 /S =0.272 ps. In Eq. (jij), the dimen- 
sionless amplitude of Morse potential e is a small parameter and is equal to 
e = D/(2Sa,dh) = 1/1400. We have also introduced the dimensionless decaying 
coefficient of Morse potential a = a^h — 7, and the dimensionless equilibrium 
distance do = cidXo/a = 0.3. 

For our analytical study, we use, as in previous numerical studies [24"! HTj. 
a Hamiltonian system composed of 30 coupled pendula. Figure [lja) is a model 
of 30 coupled pendula. Before studying its dynamics, it is instructive to look 
at a single pendulum in a Morse potential without coupling. Figure [ljb) shows 
the phase space of such single pendulum. It has two stable equilibria and a 
saddle. The black curve is the homoclinic orbit that separates two types of 
motion, namely, the oscillation near the equilibria (9 e , 0) = (±2.346, 0), and the 
flipping across the saddle (0, 0). The n-coupled pendula have similar but much 
more complicated behaviors. First, the system has two global stable equilibria, 
achieved when all the pendula have the same angular displacements 9k = 9 e 
(thus nullifying the nearest neighbors coupling) and each is positioned at the 
equilibria of a single pendulum. It also has a rank one saddle at 9k = where all 
pendula are at the upward position. For small energy, the pendula are liberating 
near one of the global stable equilibria where all angles 9k are the same and equal 
to 9 e . For large enough energy, it has been observed that a local activation can 
cause the pendula to move collectively from one energy basin to the other and 
to flip across the rank one saddle. 



2.2 Phenomenon of Structured Activations 

Previous studies [Ml [5J [HI E], mainly numerical, showed that this model ex- 
hibits an intriguing phenomenon of structured activations observed in many 
bio- molecules: while the system is robust to noise, it is sensitive to certain 
specific fine scale modes that can trigger the division. Figure [2] provides the 
numerical data for such claim. The figure shows the relationship between the 
initial amount of energy injected for various types and shapes of activation and 
the time for DNA division. 

In order to appreciate this figure and the phenomenon of structured acti- 
vations, we need to introduce the Fourier modal coordinates q which relate to 
the original system coordinates 9 through the following linear transformation 
(9 = Tq): 



2 1 

2 E 

Q = 



1 2Trka 2itka 

—7=qo + cos q a H y^qn + Sin qi 

V2 n V2 2 n i 



(3) 



where k = 1, n and a — 0, n — 1. Here, we have assume that the number 
of pendula n is even. If n is odd, we merely need to have the middle column 
corresponding to a = | removed and the column altered accordingly. 

These modal coordinates help to reveal the natural dynamics of the system 
by diagonalising the linear coupling terms and rewriting the Hamiltonian as 
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Time for Conformation Change vs. Perturbation Amplitude and Shape (mode shape) 




Figure 2: Figure shows the relationship between the initial amount of energy injected 
for various types and shapes of activation and the time for DNA division. Reprinted with 
permission from B. Eisenhower, PhD Dissertation, University of California Santa Barbara, 
(2009) ITT] . 



follows: 

n—l 71 n — 1 

H(q, P ) = {H + + e E ^(E ( 4 ) 

a=0 k = l ,3=0 

where U(9) — ( e - a ( 1 + cos6 '-' i o) — l) 2 is the Morse potential for a single pendu- 
lum. More specifically, the matrix T is nothing but the matrix of orthonormal 
eigenvectors used for the diagonalization and uj a are their corresponding eigen- 
values. By using this coordinate system, the model can be seen as a small per- 
turbation of n harmonic oscillators with frequencies ui a = yj2 — 2 cos (2ira/n), 
a = 0, n — 1. This can also be seen clearly if we write the equation of motion 
in the Lagrangian form: 

<Zo = -eM (q ,q 1 , ...,gr„_i), 

q a +u>lq a = -eM a (q , q u ...,q„_i), (5) 
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Figure 3: Figure shows a sequence of six snapshots of the evolution of 30 pendula from one 
equilibrium — e to the other equilibrium 9 e . Because the coupling is much stronger than the 
nonlincarity, the transition is collective, closely follow the mean, the thick blue line. 



where 

n n— 1 

M(q ,q 1 ,...,q n -i) = 2J u (z2 T kpqp) (6) 

fe=l /3=0 

is the Morse potential term, and Mq = dM/dqo, M a — dM/dq a are the partial 
derivatives of M with respect to qo, q a respectively. 

Observe that the coordinate of the (approximate) Oth Fourier mode, given 
as follows 

n 

8b = ^X> = V^ (7) 
fe=i 

is the average amplitude of the pendula except for a constant factor of \/n. It 
plays a special role as the collective variable, the reactive coordinate, and the 
slow variable for the system. Other (n-1) modal coordinates q a are the fine scale 
variables, the bath coordinates, and the fast variables. 

The role of approximate Oth mode (or reactive) coordinate can be seen from 
the sequence of six snapshots of the evolution of 30 pendula from one equilibrium 
—9 e to the other equilibrium 6 e of the Morse potential for a single pendulum. For 
clarity of illustration, only one pendulum is perturbed as the initial activation. 
Because the coupling is much stronger than the nonlinearity, the transition is 
collective, closely follow the mean, the thick blue line. Therefore, the average 
angle of pendula qo can be used to monitor and mark the time of DNA division. 
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Now we are ready to appreciate Figure [2j The figure shows the relation 
between the initial amount of energy injected for various Fourier modes of ac- 
tivation and the time for DNA division. The curves have been constructed by 
choosing the initial activation as a single Fourier mode. Take the magenta curve 
of the 4th mode as an example. Assuming that the initial position of the system 
is at its equilibrium point (qo, <Z( n _i)) = (— -\/n6 e , 0, 0), certain amount of 
the momentum of the 4th mode (po 5 •••iP(n-i)) — (0, 0, 0, 0,p4(0), 0, 0) is in- 
jected for the initial activation at t = and its amplitude could be modified to 
vary the amount of activation energy E a — p4(0) 2 /2. After the system evolves 
for awhile, the time of division could be determined as the time when the 
average angle go first crosses the position where go = 0. In this way, a curve 
of (E a ,td) is obtained that shows the amount of activation energy vs the time 
to DNA division for each Fourier mode. The integration time is set for 3000 
units which is approximately equal to 2/e. Notice that there are 14 such curves 
from left to right representing those from 0th mode to 13 mode respectively. 
Each curve has an asymptote at the low energy limit and will be named the 
minimum activation energy for each excited mode. White "CPs show the data 
for the random noise. From the figure, we can see that the minimum activation 
energy depends on the way this energy is injected into the system. While the 
system is robust to noise, it is sensitive to certain specific fine scale modes that 
can trigger the division. There may also exist a relation between the minimum 
activation energy of each Fourier mode and its modal frequency. In this paper, 
we would like to develop the analytical methods to reveal the activation mecha- 
nism and to compute the minimum activation energy for each mode. Moreover, 
we want to develop the techniques for controlling the real DNA division via low 
intensity electro-magnetic fields and to reveal how enzymes initiate the DNA 
opening dynamics. 



3 Analytical Study of Structured Activations 

Let us start our analytical study of the phenomenon of structured activations. 

3.1 Nearly : 1 Resonance and Partial Averaging 

Recall that the equations of motion of our system are given by Eq. |5]) which can 
be seen as a nonlinear perturbation of n harmonic oscillators with frequencies 
uj a = y^2 — 2 cos (2-7ra/n), a = 0, 1, n— 1. Note that besides luq = 0, uj a varies 
from 0.2091 to 2. Hence, the reactive mode forms a nearly : 1 resonance with 
any other mode, each of which has an 0(1) frequency: 

muiQ + Qui a = 0, with m = 1. (8) 

While such a relationship may not be seen as a resonance in the classical "engi- 
neering sense" , it certainly is in the mathematical sense. This fact leads to small 
denominators and coupling terms in the corresponding averaged equations or 
normal form. Other modal frequencies, from 1st to 14th mode (and from 15th 
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to 29th mode), are not rationally commensurate and do not have significant 
time scale separation. We do not expect strong resonance among them. Hence, 
we believe that nearly : 1 resonance should be the main focus of our study. 

Nayfeh et al. [23 IMj, Haller Q2], Feng and Liew QJ], Tuwankotta and 
Verhulst [35], Langford and Zhan j^Hj and Broer et al. [5] have studied such 
a degenerate resonance. Both Langford and Zhan and Broer et al. used the 
method of normal form. Nayfeh et al., Haller, and Feng and Liew applied 
a modified averaging method directly to a two mode truncation of a simple 
mechanical system with parametric or external excitation. Tukwankotta and 
Verhulst applied a similar method to the study of nonlinear wave equations. 

While Eisenhower and Mezic did not mentioned the term, nearly : 1 res- 
onance, in their papers, they did apply Arnold's method of partial averaging 
[I] to a truncated Hamiltonian in the study of a chain of Duffing oscillators. 
In [12] . the reason that Eisenhower and Mezic gave for studying the Duffing 
case was that "The exponential form of the Morse potential makes analytical 
progress difficult,..." Since we have not seen this kind of simplification in the 
DNA literature, we will, in this study keep the Morse potential (and the added 
difficulty) . 

3.2 Two Mode Truncation Is Adequate 

In order to carry out the analytical work for such a high dimensional system, 
certain reduced model is needed. The success of Eisenhower and Mezic [12], Du 
Toit et al. 0, Nayfeh et al. [22 EH], and other researchers [22 [181 [35] shows 
that the method of modal truncation is a viable method if it is used with care. 
Our understanding of the nearly : 1 resonance in our DNA model and the 
results of our extensive numerical simulation have convinced us that two mode 
truncation is adequate for our analytical study of the phenomenon of structured 
activations. 

Figure [4] shows the projections of a sample trajectory on the phase spaces 
of the first 15 modes. The activation has been chosen to be a single 6th mode. 
From the numerical simulation which use qo(0) = \fnQ & and P6(0) = \J1E as 
the initial condition (where E ~ 1.221 is the amount of energy injected into the 
6th mode), we observe that (i) the energy transfer takes place over-whelmingly 
from the excited mode to the reactive mode, inducing the division, (ii) only an 
extremely small amount of energy transfers from the excited mode to one or two 
other modes via near resonance. Hence, we expect that the two-mode or three- 
mode truncations i.e., a truncated system that includes the reactive mode, the 
excited mode, and perhaps a mildly affected mode, should provide an adequate 
reduced model for our analytical study on the nearly : 1 resonance, the energy 
transfer, and the activation mechanism. 

What follow are the equations for any two-mode truncation (the reactive 
mode and the 7th mode) 



eM^(qo,q y ) 
eM^(q ,q y ) 



(9) 
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Figure 4: The fi gure shows the projections of a sample trajectory on the phase spaces of 
the first 15 modes. E = 1.221 is the amount of initial kinetic energy injected into the 6th 
mode. 



where 

n 

M( 2 )(g ,g 7 ) = ^C/( Y, T ^P) ( 10 ) 

k=l /3={0, 7 } 

is the reduced Morse potential term, and M = dM ( ^ /dq , M^ 2) = dM^/dq-y 
are its partial derivatives with respect to qo, q 7 respectively. Note: for notational 
simplicity, M(q ,q 7 ) will be used later for M^ 2 \q ,q 1 ) whenever there is no 
ambiguity. 

3.2.1 Error Estimates for the Two Mode Truncation. 

Here we show that the solution trajectories of our two mode truncation Eq. (|9| 
are within 0(e) of the solution trajectories of the original full system described 
in Eq. Q for at least 0(1) times. 

The proof is an extension of the one in Du Toit et al. [5] , where a general 
system that has the same form as Eq. ^ was studied. In that paper, the 
authors (i) proposed a general technique for obtaining an 1 ^ degree of freedom 
reduced system and (ii) were able to provide a rigorous error estimate for their 
procedure. First, they introduced an approximation by replacing Eq. [5jb) with 
the analytical solutions of the unperturbed linear system, as defined by 

Qa = A a cosui a t H smu! a t, a = l,...,n— 1 (11) 

UJ a 
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where A a ,B a are the initial conditions. Then, they obtain their reduced equa- 
tion 

Q = -eM (Qo,Qi,....Q n -i) (12) 

via simple substitutions. Hence, the information contained in the other modes 
persists in the reduced equation of the reactive mode via the initial conditions. 
The error estimate that they provided claims that the solution trajectories of 



Eq. (12) are within 0(e) of the solution trajectories of the original full system 
described in Eq. ^ for at least 0(1) times. This result is shown by applying 
a standard error analysis technique: substituting a formal expansion of the 
solutions, using the Lipschitz continuity of Mo, and then applying the Gronwall 
lemma as is done in the proof of Theorem 9.1 in [37J (see also [33] for a more 
detailed analysis on why removing small nonlinear perturbation to harmonic 
oscillations incurs small error). 

Notice that the equations of our two mode truncation Eq. ^ can be ob- 
tained by simply setting all the initial conditions A a , B a for Eq. ^ equal to 
zero except if a = 7. Hence, the solution trajectories of our two mode trun- 
cation Eq. |9]) should also be within 0(e) of the solution trajectories of the 
original full system described in Eq. ^ for times 0(1). 

Moreover, this error estimate is valid for an arbitrary forcing function. But 
for our DNA model, the exponential decay of the Morse potential with the 
distance implies that when the pendula escapes the immediate vicinity of the 
opposing pendula, the Morse potential and its consequent perturbation are ef- 
fectively zero. Hence, our error bound is very loose because it utilizes only the 
fact that e is small, but not the specific feature of our DNA model, which is 
that the forcing term is also small in a large region of phase space. This is why 
we numerically observed that the two mode truncation remains accurate over a 
timescale much longer than 0(1). 

Figure [5] show the amount of activation energy vs the time of DNA division 
for the 6th mode. The blue "o"s are the data points for the full equations of 30 
modes; the red "*"s are correspondent data points for a two mode truncation 
(the reactive mode and the 6th mode). Both have almost the same minimum 
activation energy 1.205. Even their time of division differs very little if the 
actuation energy is slightly larger than the minimum activation energy (away 
from the asymptote). The above observation also holds true for all the other 
modes. 

Now we are ready to do the truncation and apply the method of partial 
averaging to obtain the averaged equations for the reduced models of this chain 
of coupled Morse oscillators, and use them to study the activation mechanism 
and compute the minimum activation energy for each mode. 

3.3 Averaged Reduced Hamiltonian 

Since in our case, the partial averaging of Euler-Lagrangian equations is equiv- 
alent to the partial averaging of its Hamiltonian and our main concern is on 
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Figure 5: Figure shows the amount of activation energy vs the time of DNA division for 
the 6th mode. The blue "o"s are the data points for the full equations of 30 modes; the red 
"*"s are correspondent data points for the two mode truncation. Both have almost the same 
minimum activation energy 1.205. 



energy transfer, we preferred to use the Hamiltonian formulation in this study. 
The equivalency was shown in ^jA] of the Appendix. 

Recall that the Hamiltonian of the two mode reduced model is given by 

n 

Hl(q,p)= ]T (H + H<&)+*Y, U ( E T WP)- ( 13 ) 

"={0,7} fe=l /3={0,7> 

The averaged reduced Hamiltonian can be obtained as follows. First, we ap- 
proximate (via Taylor expansion) the Morse potential U, which involves the 
exponential function, by a polynomial of degree 26 at 9 = 

n 26 

Q={0,7> k=lj=0 £={0,7} 

where a,j,j — 0, ...,26 are the coefficients of the expansion. As pointed out in 
[12], the Morse potential is indeed a difficult function to work with. Since the 
geometry of the Morse potential requires a polynomial of degree at least 26 for 
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an accurate approximation, we decide to use such a high degree expansion for 
our computation of the averaged reduced Hamiltonian. 

Then, we use the angle-action coordinates defined as follows 

<7 7 = J 2/ 7 /cl> 7 sin 7 , p 1 = ^J2I 1 uj 1 cos (/> 7 . (15) 

and rewrite the reduced Hamiltonian as iTJ (<Zo,??o, Ij,^)- Notice that, besides 
qo,Po, the action I 1 is also a slow variable. Hence, the averaged reduced Hamil- 
tonian can be obtained by averaging the only fast variable 7 : 

/■27T 

^2(90,^0,^7) = ^/ #2(?0,P0,-f T ,</M# 7 - ( 16 ) 
JO 

After renaming variables, the averaged reduced Hamiltonian is given by 

H 2 = ^y 2 +ojI + e (na + £ c 2k (I)xA (17) 

V k=0 ) 

where x = q~o,y = po are the Cartesian coordinates of the reactive mode; I, u 
are the action and the constant frequency of the other mode; enao = 0.0214 is 
the energy value at the saddle; c%k(I) are polynomials in /. For example, Cq(I) 
is a 13 degree polynomial in / given by 

13 

c (I) = J2bjP (18) 
3=1 

where bj are numerical coefficients. Also, for the simplicity of notation, the 
letter 7 is dropped from the averaged reduced Hamiltonian ■ 

3.3.1 Error Estimates for the Averaging. 

By applying the standard averaging theory for differential equations (see for 
instance |30j : we again recall that averaging the Hamiltonian is equivalent to 
averaging the equations), we can conclude that the solution trajectories of the 
averaged reduced equations are within 0(y/e) of the solution trajectories of the 
reduced equations for at least 0(1/ y/e) times. 

3.4 Error Estimates of the Entire Treatment. 

There are three approximations that we made in our treatment: Taylor ap- 
proximations of the potential, truncation into two modes, and averaging. We 
found that the truncation and averaging respectively induces 0(e) and 0(^/e) 
errors for at least 0(1) and 0(1/ ^/l) times. In addition, as long as the solution 
remains bounded, Taylor approximations result in a small o(e) error in the non- 
linear forces, which again by Gronwall only induces 0(e) error in the solution 
till at least 0(1) times. Put together, our entire treatment induces at most an 
0(y/l) error till at least 0(1) time. 
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Fi gure 6i (a) The contour plots for the averaged reduced Hamiltonian when 7 = 5. This 
phase space (x, y) has a separatrix which separates two types of motion: liberation (also 
referred to as "breathing" in [2]) and flipping, (b) Three such phase spaces stacked up in 
increasing I in the (x,y,I) space. Notice that the phase space and the separatrix "shrinks" 
towards the saddle as I increases. Together, these separatrices form a homoclinic manifold 
and can be used in studying the minimum activation energy for each excited mode. 



Moreover, this error estimate is valid for an arbitrary forcing function. But 
for our DNA model, the exponential decay of the Morse potential with the 
distance implies that when the pendula escapes the immediate vicinity of the 
opposing pendula, the Morse potential and its consequent perturbation are ef- 
fectively zero. Hence, our error bound is very loose because it utilizes only the 
fact that e is small, but not the specific feature of our DNA model, which is that 
the forcing term is also small in a large region of phase space. This is why we 
numerically observed that our approximation remains accurate over a timescale 
much longer than 0(1). 

3.5 Phase Space of Averaged Reduced Equations 

Given the averaged reduced Hamiltonian, we obtain the averaged reduced Hamil- 
tonian equations as follow 

x = y, 

y = -e(f^2kc 2k (I)x 2k -Ax; <j> = uj + e ( £ ^z 2 ' j . ( I !> i 

Notice that the action I of the excited mode is a constant of motion. Hence, 
the averaged equations of the reactive mode, i.e., the x,y equations, and their 
phase space structures are parametrized by i\ 




Therefore, the averaged reduced Hamiltonian, Eq. (17), can be used to study 



the phase spaces of the reactive mode of the averaged reduced equations that 



14 



are parametrized by /. Figure |6[a) shows the contour plots of the averaged 
reduced Hamiltonian and the phase space for the reactive mode of the averaged 
reduced equations, when 1 = 5. Notice that this phase space has a separatrix 
which separates two types of motion: liberation and flipping. Figure [6jb) shows 
three of such phase spaces stacked up in increasing I (action of excited mode). 
In the (x, y, I) space, the vertical axis can also be seen as the axis of increasing 
activation energy by scaling with w, E act = loI. Notice that the phase space 
and the separatrix "shrinks" towards the saddle as E act increases. Together, 
these separatrices form a homoclinic manifold and can be used in studying the 
minimum activation energy for each excited mode. 



3.6 Analytical Study of Minimum Actuation Energy 

Recall that for large enough energy in the excited mode, the energy transfered 
to the reactive mode will surpass those at the saddle and the system will be 
driven from its initial equilibrium across the separatrix, causing the flipping. 
See Figure |4| But for the averaged reduced system, this process manifests 
itself via the changes in the phase space and the separatrix of the reactive 
mode parametrized by /. See Figure j6^b) for illustration. For activation energy 
slightly larger than the minimum activation energy, E m i n (= toI m ), the initial 
stable equilibrium at (x e , 0), now parametrized by I^ in , will cross the separatrix 
parametrized by J+, move into the flipping region of phase space parametrized 
by and induce the DNA division. Hence, E min can be found by the condition 
that the point (x e ,0,I m ) is on the separatrix that passes through (0,0, I m ). 
Since the separatrix is the energy curve (x, y) defined by 

H 2 (x,y,I m ) = H 2 (0,0,I m ) (20) 

where H 2 is the averaged reduced Hamiltonian, E m i n can be found by solving 
the following equation for I m 

H 2 (x e ,0,I m )=H 2 (0,0,I m ) (21) 

where (0, 0, I m ) is a saddle at I m and {x e ,0) is the initial equilibrium point at 
1 = 0. 



As pointed out earlier, H 2 of Eq. (17) has been computed using a Taylor 



expansion of the Morse potential U at 9 = 0. Strictly speaking, it may be better 
to denote it as H®- Even though we have already tried to find an excellent 
Taylor expansion that can take care of the approximation of U at both the 
neighborhood near the saddle 9 = and the basins around the stable equilibria 
9 = ±9 e , the accuracy of U and iJ^ at the latter is still not as good as those 
at the former. For example, while the actual stable equilibria should be at 
(±Vn0 e ,O) = (±12.85,0), the approximated ones will be at (±12.59,0) if H% 
is used. Since any significant error in the approximation of the Morse potential 
U at the stable equilibria ±0 e will cause a large error in the estimation of the 
minimum activation energy, another averaged Hamiltonian iff has also been 
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Frequency-Minimum Actuation Energy Curve (cyan o: numerical; red *: analytical) Frequecy-Mininum Energy curve (Full Eq vs 2Mode Truncation) 
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Fi gure 7'. (a) Analytical estimation of minimum activation energies vs numerical simulation 
obtained from two mode reduced models. The data for magenta "*"s are from analytical 
computation. The data for cyan "o"s are from numerical simulation obtained from the two- 
mode truncations. They match very well, (b) Numerical simulation of two mode reduced 
models vs numerical simulation obtained from the full model. The blue "o"s are the data 
from numerical simulation of the full equations (30 modes). The red "*" are from simulation 
of the two-mode truncations. 



computed using an expansion of the Morse potential U at one of the stable 
equilibria 9 = —6 e : 

H! = \y 2 + ljI + e (y^c'jWix + Vn9 e )^j . (22) 

Here, the letter e stands for the word "equilibria". Now, The minimum activa- 
tion energy can be analytically estimated by solving the following equation 

fff(-Vn0 e ,O,/ m ) = ff£(0,0,/ m ). (23) 

Straightforward substitution will give us 

ul m + ec' (I m ) = U)I m + e(na + c (I m )). (24) 

After simplification, the minimum activation energy can be analytically esti- 
mated by solving 

c' {I m ) - c (I m ) - na = 0. (25) 

For the two-mode truncation of the reactive mode and the 6th mode, I m /ujQ — 
0.8539. Therefore, E min — uj 6 x 0.8539 x cj 6 = 1.1801, as compared with 1.205 
from numerical simulation. 

Figure[7]Ja) compares the analytical estimation of minimum activation energy 
with the numerical simulation obtained from two-mode reduced model. The 
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data for cyan "o"s are from the numerical simulation obtained from the two- 
mode truncation. For example, the 6th data point (for the 6th mode) with 
uiq = 1.17 needs minimum activation energy = 1.205. The data for magenta 
"*"s are from analytical computation. For example, the 6th data point with 
loq = 1.17 has E^ in = 1.1801 Clearly, the data from numerical simulation and 
the data from analytical computation match very well. Moreover, E m i n (except 
the 15th mode) can be approximated as a parabola 

E mm = 0.8539 x uj 2 . (26) 

Here, we would like to remark that the use of two different Taylor expansions 
of the Morse potential function U, one at 6 = and another at 9 = 6 e , not only 
allows us to deal with the difficulty posed by the exponential form of the Morse 
potential function but also improves the accuracy of our analytical estimation 
of the minimum activation energy for each excited mode. 

Figure [7](b) shows the comparison between numerical simulation obtained 
from the two mode truncations with the numerical simulation obtained from 
the full model. The blue "o"s are the data from numerical simulation of the 
full equations (30 modes). For example, the 6th data point (for the 6th mode) 
with loq = 1.17 needs minimum activation enesrgy = 1.205. The red "*" 
are from simulations of the two-mode truncations Em = 1.205. Besides the 
14th mode (with error less than 5%), all other data values of blue "o"s and red 
"*"s have differences less than 1%. Therefore, for the study of the minimum 
actuation energy, the analytical estimation provides accurate prediction for the 
full system. . 



3.7 Summary of Analytical Results on Structured Activa- 
tions 

By applying the method of partial averaging, we have obtained the averaged 



reduced equations Eq. ( 19 1 for a chain of coupled Morse oscillators, (i) These 
equations reveal the coupling between the action and energy of the excited mode 
and the dynamics of the reactive mode, as well as the phase space structure of the 
activation mechanism, (ii) They allow us to estimate analytically the minimum 
activation energy for each excited mode and discover a relation between the 
frequency of the excited mode and their corresponding minimum activation 
energy, (iii) These estimates match very well with the numerical simulations 
obtained from the reduced and the full model. The results show that the nearly 
: 1 internal resonance is responsible for the structured activation of our DNA 
model. 



3.8 Remark on Pitchfork Bifurcation 

We observe that the reactive mode of the averaged reduced model has a pitchfork 
bifurcation at lb ps 38. This can be seen clearly by studying the following 
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Figure 8: Figure (a) shows that the reactive mode of the averaged reduced system has a 
pitchfork bifurcation at the parameter value / 38. Figure (b) shows a graph of 02(1) which 
has a simple zero at I « 38. 



averaged equations of the reactive mode parametrized by /: 

x = y, y=-e ^ 2kc 2k {I)x 2k ~ 2 ^j x; (27) 

Notice that the equilibria of this system are given by (0,0) and (x e (I),0) 
where x e (I) are the solutions of 

13 

Y,2kc 2k (I)x 2k - 2 =0. (28) 
fe=i 

See the blue curve in Figure ^a,). Clearly, (x e {I), 0) are stable equilibria, except 
at (0, 0). In fact, they are surrounded by liberation contours as shown in Figure 
[6j As for the equilibria (0, 0) for each /, their stability are determined by the 
sign of the coefficient of the linear term in x, namely, — 2ec2(/). For C2{I) < 0, 
they are saddles. For 02(1) > 0, they are stable equilibria. Figure [8](b) shows 
a graph of 02(1) which has a simple zero at lb ~ 38. Hence, the system has a 
pitchfork bifurcation as shown in Figure^ a). As will be shown in the following 
section, this bifurcation will play an important role in our study on the control 
of this DNA model via parametric resonance. 

4 Control via Parametric Resonance 

Building on our understanding of its internal dynamics, we want to control the 
division of this DNA model via parametric excitation, that is in resonance with 
its internal trigger modes. This effort is guided by two observations and two 
conjectures: The two observations are (i) the averaged reduced model has a 
pitchfork bifurcation at lb ~ 38 — the physical interpretation of this is that 
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if enough energy is injected into the trigger mode with a value larger than its 
bifurcation value, the chain of pendula will remain near its open state; (ii) para- 
metric resonance is an efficient way for energy transfer from an external source. 
Hence, our two conjectures are (i) low intensity electro-magnetic fields can be 
used with the mechanism of parametric resonance to pump energy into the trig- 
ger modes and keep the real DNA chain near the opening state for replication 
and transcription. Here we note that electric field-induced molecular vibrations 
have already led to the development of a noninvasive cell transfection proto- 
col that enables foreign DNA molecules to cross cell membranes and penetrate 
into the cytoplasm by eliciting vigorous vibration between molecules and cells 
|44l 133] ; (ii) enzymes may use similar method to initiate the opening dynamics 
for DNA replication and transcription. For discussions and evidence of impor- 
tant relations between low-frequency vibrations in DNA molecules (and other 
bio-macromolecules) and significant biological functions, we refer to [4], [5], [3], 

m- m- cm m- m 

Our three main results in this section are as follows: (i) We identify excita- 
tion parameters (as a function of friction) and a class of trajectories that show 
how parametric resonance can be used to drive the averaged reduced model 
from its almost equilibrium state to its open state, (ii) By identifying the av- 
eraged reduced effective Hamiltonian (after parametric excitation), we uncover 
the global phase space structure and analyze the characteristic trajectories men- 
tioned above, (iii) We extend the results for the averaged reduced model with 
parametric excitation and friction to the reduced and to the full model with 
parametric resonance and friction. The findings uncover a method for control- 
ling DNA division via parametric resonance. Although the issues of inhomo- 
geneity, helicity, and environmental effects (such as noises) are ignored for now, 
they will be taken into consideration in our future work. 

4.1 Equations of Motion with Parametric Excitation 

The equations of motion of our full DNA model with parametric excitation and 
frictions can be written as follows: 

h - (0 k +i - 29 k +0 k -i)- eU'(8 k ) = e8 k f cos Qt - efi9 k (29) 

where k = 9q = 9 n , U is the Morse potential function and U' is its 

derivative. Notice that the left hand side is the original equations of motions, 
without parametric excitation or frictions. See Eq. [2] Moreover, /, SI are the 
amplitude and the frequency of the parametric excitation, respectively; SI is in 
a nearly 1:2 parametric resonance with a chosen internal trigger mode w 7 ; [i is 
the frictional coefficient. 

Note the friction terms represents energy loss caused by interaction with 
surrounding molecules. Noise terms (ignored here) would represent energy gain 
caused by (thermal) interaction with surrounding molecules [25 . The paramet- 
ric excitation terms represent interactions with electro-magnetic fields, we note 
that DNA vibrations induced by electric-fields or microwave absorbtion are an 
experimental reality [44] [33] . 
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As pointed out in Ref. [H], the bases of DNA have dipole moments [7], 
which could couple with an external electromagnetic field. In principle, this 
coupling can induce a periodic wave-like forcing over DNA bases. Since an 
electromagnetic field has a polarization, it is reasonable to assume that the 
coupling between an electromagnetic field and a DNA base is dependent on 
the orientation of the base, Ok- This motivates us to introduce the parametric 
excitation term in our model equation (the first term on the right-hand side of 



Eq. {29) above). 

After using the Fourier modal coordinates, the equations of motion for a 
corresponding two mode reduced model are given by 

qo = Po Po = -eM - efipo + efqo cos fit. 

q~ f = p 7 p-y = —v^qry — eM 7 — efip 7 + efq 7 cos fit (30) 

where qo,Po and q-y,P*Y are the coordinates for the reactive mode and the 7- 
mode, respectively; M is the reduced Morse potential term and Mq,M 7 are its 



partial derivatives with respect to qo,q^ respectively. See Eq. (|9|) and Eq. (10 1 
for detail. 



After applying the method of partial averaging to Eq. ( 30 ) using the angle- 
action variables 



In 



= y/4I/flsm(%t + p) , Vl = V70cos(ft + /3) , (31) 

renaming variables, and setting fl = 2uj, we obtain the averaged reduced equa- 
tion of motion 

x = y, ft = e(M/ - a/2u + / cos 2/3/4w), 

y = -e(M x + »y), t = el(f sin 20/2w-p), (32) 

where J, (3 are the action and the phase of the chosen internal trigger mode; a 
is a detuning parameter, defined by fi 2 /4 = cj 2 + ea; M is the averaged Morse 
term, and M x ,Mj are its partial derivatives with respect to x,I respectively. 



See Eq. ( 19 ) for comparison. 

The derivation is similar to what has been done in [551 H] an d is parallel to 
the method used in fAlof the Appendix with only minor modification. We first 



rewrite Eq. ( 30 ) in the following first order form: 

qo = VePo Q-y = P~ t 

Po = -VtRa Pj = -u^q-y - tR-f. (33) 

where 

Ro = M + ppo - fqo cos fit 

R~! = M 7 + np-y - fq 1 cos fit. (34) 



By using the angle-action variables defined by Eq. (31), we can transform Eq 



( 33 ) into the Lagrangian Standard Form 

qo = Vepo P= e{R 1 /^/lfi)si\\(j) 

po = -Vei? I = — cRy yj il /fl cos <f> (35) 
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where <f> = (jt + 0) and 

fi 7 = i? 7 - a V \l JCl sin 4>. (36) 

Hence, we can apply the standard averaging theory (by averaging t from to 
47r/f2) and obtain the averaged reduced equations 

e f 2n 

x = y/ey $ = — (Av/\/M)sin 

2tt Jo 



y 



7/ 

2?r , /•2ir 



x ' R d^ i = -~l (RyV^I/Si) cos <f>d<j> (37) 



2tt ./o 2?r Jo 



where x = (?0j V = Po- Then, after carrying the averaging (that is similar to jjA|, 
renaming variables, and setting Q — 2w, we obtain Eq. ( p2| ). 

Moreover, the averaged reduced system (without friction) has an effective 
Hamiltonian 

H PB ^ + ta -,f J + ^ ,38) 

that can provide insights on the global phase space structure of this averaged 
reduced model. Notice that since angle-action variables are used in our deriva- 
tion, the Hamiltonian that we have obtained is canonical. See reference [IB] for 
comparison. 

4.2 Merging Local Bifurcation Analysis with Global Ge- 
ometry of the Effective Hamiltonian 

Detailed local bifurcation analysis of the averaged reduced system can be used 
to reveal the ranges of a and other parameters, /, /z, where the desired dynamics 
may be available [26l [28]. For example, for / = 2.5, /i = 0.5/ai where u = ujq — 
1.17, the frequency response curves for the averaged reduced system can be 
depicted as in Figure [10] These frequency response curves can be obtained by 



studying the fixed points of Eq. ( 32 ) and their stability that are parametrized 
by the detuning parameter a . Notice first that for all fixed points, y = y e = 0. 
Moreover, there are two types of fixed points, / = or / ^ 0. 

4.2.1 Fixed Points with 1 = 0. 

There exist two sub-cases: (i) x = and the fixed point is (x, y, q,p) = (0, 0, 0, 0), 
(ii) x = x e = ±12.59 and the fixed points are (x,y,q,p) = (x e , 0,0,0). These 
two types of solutions can be obtained by solving the following equation 

M x = 2kc 2k (I)x 2k -^j x = (39) 

with 1 = 0. Notice that in our study of the fixed points and their stability with 
7 = 0, the following Hamiltonian polar coordinates has been used: 

q = \Jll /u) sin 13, P = V2I oj cos (3 (40) 
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That is why (q,p) = (0,0) when 1 = 0. 

For case (i), the eigenvalues of the fixed point are given by 

A= - £/i± ^ )2 - 8£b , -f ±^VP-^-2by (41) 

where b = C2(0) = —0.47 x 10~ 4 . For example, as a decreases from 3 to 
—2, this fixed point goes through two bifurcations, one at 1.1455 and an- 
other at —1.1457. In the interval (—1.1457,1.1455), it is a rank two saddle 
(saddle x saddle), marked with red dashed lines. In the intervals (— oo, —1.1457) 
and (1.1455, oo), it is a saddlexstable foci, marked with magenta dashed and 
dotted lines. See Figure [l0|a) (b) . 

For case (ii), the eigenvalues of the fixed points are given by 



X= -*» ± ^f - 4£ M ^, -^±yp^ { a-2uM lY (42) 

where M xx is the second partial derivative of the averaged Morse potential term 
M with respect to x and Mj is the first partial derivative of M with respect to I. 
In the interval (—2, 3), these two fixed points (±12.59, 0, 0, 0) are stable (stable 
focixstable foci), marked with blue lines. See Figure [l0^c)(d). Even though 
they do go through bifurcations at 55.6726 and 57.9636, we do not include their 
analysis in this paper because they are far outside of the region where we have 
found the desired dynamics. 

4.2.2 Fixed Points with 1^0 

For these fixed points, /3 = e = | arcsin y / 2w^//. Again, there exist two 
sub-cases, (iii) x=0 and the fixed points are (x,y,I,0) = (0, 0,/ e ,/3 e ) where 
le = lei ') i s obtained by solving the following equation for I e 



M > 2. 



l -[a±\y/P-^)^j (43) 



with x — 0. See Figures |9ja) (b) for the curves of these fixed points, (iv) x ^ 
and the fixed points are (x, y, I, 0) = [x* , 0, I*, f3 e ) where I* = I*[a),x* = x*(a) 



are obtained by solving the two nonlinear equations Eq. ( 39 ) and Eq. ( 43 1 
simultaneously. See Figures [9](c) (d) for the curves of these fixed points. 
For case (iii), the eigenvalues of the fixed points are given by 



A= -g^F4g ; -f ±^^) 2 + 4.//M // cos2/ 3 , (44) 

Numerical computation shows that (1) the fixed points of the upper left branch 
of I e = I e {u) in Figure ^ a) are stable (stable focixstable foci) for —1.4825 < 
a < 1; (2) the fixed points of lower left branch are saddle x saddle; (3) the 
fixed points of right branch are saddlexstable foci. Note: for simplicity of 
presentation, only part of the x — solution in Figure |9jb) that corresponds to 
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Fixed Points with J=J and x=0 



Fixed Points with J=J and x=0 




Figure 9: Figures (a)(b) and Figures (c)(d) show the curves of fixed points for I e = 
I e (a),x = and I* = I*(cr),x* = x*(cr) respectively. Figure (a): (1) the fixed points of 
the upper left branch of I e = I e {cr) are stable (stable foci X stable foci) for —1.4825 < a < 1, 
marked with the blue line ; (2) the fixed points of lower left branch are saddle X saddle, marked 
with the red dashed line; (3) the fixed points of right branch are saddlexstable foci, marked 
with the magenta dashed and dotted line. Figure [9jb) : for simplicity of presentation, only 
part of the x = solution that corresponds to upper left branch in Figure (a) has been marked 
with the blue line (stable focix stable foci). Figures (c)(d): (1) the fixed points of the left 
branches are saddlexstable foci, marked with the magenta dashed and dotted lines; (2) the 
fixed points of the right branches are stable (stable focix stable foci), marked with the blue 
lines. 



upper left branch in Figure |9]ja) has been marked with blue line (stable focix 
stable foci). 

For case (iv) , the characteristic equation of the eigenvalues of the fixed points 
is given by 

A 4 + tiA 3 +t 2 A 2 +t 3 A + t4 = (45) 

where 

T\ = 2e/i 

r 2 = e 2 fi 2 + eM xx - e 2 fIM n cos 2/3/w 
r 3 = e 2 fiM xx - e 3 fifIM n cos 2(3 /uj 
r 4 = e 3 fI{M Ix - M xx M n ) cos 20/u. 
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Fi gure 10i Frequency response curves. Figure shows how the fixed points and their 
stability change as the detuning parameter cr is varied. Since for all the fixed points, y = 
and /3 = /3 e , only (a) I vs a and (b) x vs <r are plotted. For example, in the neighborhood of 
a = (—0.5, 0.5), the system has four types of fixed points for each <r: (1) the red dashed curves 
where I = 0, x = denote saddle X saddle; (2) the solid blue curves where I = 0, x = x e = 
±12.59 denote stable focixstable foci (only positive x e is drawn); (3) the solid blue curves 
where I = I e ,x = denote stable focixstable foci; (4) the magenta dashed dotted curves 
where I = I*, x = x* denote saddlex stable foci. 

Numerical computation shows that (1) the fixed points of the left branches of 
Figures |9jc) (d) are saddlexstable foci; (2) the fixed points of the right branches 
are stable (stable focixstable foci). 

Figure [TO] is the combined result of the case studies. The frequency response 
curves in this figure show how the fixed points and their stability change as the 
detuning parameter a (and hence the frequency Q) of parametric excitation is 
varied. By analyzing the relationship between all these curves, we can see that 
the averaged reduced model may have the desired dynamics in the neighborhood 
of a = 0. For example, at a = 0, the phase space has four types of fixed points 
as follows: 

1. (x, y 7 q,p) — (0, 0, 0, 0) is a rank two saddle. 

2. (x,y,q,p) = (x e , 0,0,0) where x e = ±12.59 are stable foci. Notice that 
these fixed points locate at the bottom of the potential well of the averaged 
reduced system. 

3. (x,y,I,/3) = (0, 0, I e ,Pe) where I e = 45.74 x w 6 ,^ e = \ arcsin y/2ijj 6 fi/f 
are stable foci. Notice that the (x, y) coordinates of these fixed points 
mark the DNA division. 

4. (a;, y, 1, 0) = (x*, 0, /*, j3 e ) where x* = ±8.7478, I* = 5.312 x lo 6 are stable 
foci x saddle. If the frictional coefficient /i is small. These fixed points 
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Fi gure 11; (a) Contour plots of the effective potential energy in the (x,I) space, (b) An 
example of a trajectory that shows how the parametric resonance drives the averaged reduced 
system from its almost equilibrium state to its open state in the (x,y,I) space. 



are essentially rank one saddles. 

According to the theory of tube dynamics [TH1 |T51 [TB] , these modified rank one 
saddles may provide a low energy pathway from the neighborhood of the bottom 
of the well to the region which marks the DNA division. 



4.2.3 Global Geometry of the Effective Hamiltonian 

While the local bifurcation analysis does provide many basic ingredients for our 
study, it does not by itself give a clear and global picture of the dynamics of 
the averaged reduced system. Hence, the effective Hamiltonian Hpf>(x,y, I, f3) 
is needed to fill in this gap. Notice that in the Hamiltonian polar coordinates 



Eq. (40), j3 = 7r/2 corresponds to the case where p = 0. Therefore, 



H PR {x,Q,I,Tt/2) = eM-e^I-ef^ (46) 
provides an effective potential for the averaged reduced system. 



Figure 11 a) shows the energy contours of this effective potential when a = 0. 
It provides us with the insights for the global phase space structure of the 
averaged reduced system. From the figure, we can see clearly how the four types 
of fixed points obtained previously fit together within the global geometry of 
the effective Hamiltonian. Moreover, 

• The parametric excitation represented by the parameters /, a has turned 
(0, / e ), which marks the DNA division, into a sink, 

• it also creates two low barrier rank one saddles that are close to the two 
DNA equilibrium states (±12.59,0). 



25 



Hence, the addition of parametric excitation to the averaged reduce system 
should allow certain trajectories with a little energy in the trigger mode to 
move from an almost equilibrium state, over the saddle, navigate down the 
energy contours, and reach the sink. All these insights drawn from the local 
bifurcation analysis and the effective Hamiltonian enable us to generate a class of 
trajectories that show how the parametric resonance drive the averaged reduced 
system from its almost equilibrium state to its open state. Figure [Ti|b) shows 
an example of this kind of trajectories in the (x, y 7 1) space. 



4.3 Parametric Resonance Drives DNA to Division 

The data for this trajectory are given as follow. For the system parameter, 
we have / = 2.5, /i = 0.5/ujQ,a = 0. For the initial condition, we have x = 
12.59,yo = 0, Iq = 0.54015 x u; 6 ,/3o = 0. The integration is done using the 



averaged reduced equations, Eq. (32). Notice that this trajectory in Figure 
|ll[ b) starts at the equilibrium position of the reactive mode but with certain 
small amount of energy in the trigger mode (6th mode). Without the parametric 
excitation, the system will liberate near the equilibrium state if there is no 
friction or die down if the friction exists. However, if the parametric excitation 
is turned on at t = with the data provided above, the 1 : 2 parametric 
resonance will inject energy into the trigger mode, increase the value of /, make 
the trajectory to reach the region that marks the DNA division (x = 0) but 
with large energies in the trigger mode. See Figure [14] below for a physical 
interpretation of this class of trajectories when it is near the DNA open state. 

Here, we would like to make a remark on the amount of initial energy in 
the trigger mode. Note that if we increase the amplitude / of the parametric 



excitation, the magenta dashed dotted curve of Figure 10 'a) will shift downward 



Similarly, the rank one saddles in Figure 11 ^a) will also shift downward. These 
mean that the amount of initial energy needed in the trigger mode, namely, the 
value of /, can be lower for large /. Numerical simulations of the full system 
confirm this observation. 



4.4 



Extend the Results to the Reduced and the Full Mod- 
els 



Figure 12 shows two corresponding trajectories, one for the reduced model and 
another for the full model. The trajectory in Figure [T2][ a) is generated with the 
following data. For the system parameters, we set / = 2.5, fj, = 0.5/uie, a = as 
before. For the initial condition, we set x = qo(0) = 12.85, y = Po{0) = 0, J = 
0.8680 x W6,/?o = and integrate the trajectory using the reduced equations, 
Eq. (30) where 



i(0) = V2/oMsinfl,, 



p 6 (0) = \/2Iouj 6 cos/3 - 



(47) 



As for the trajectory in Figure 12 Tb), it is generated with the following 
data. For the system parameters, we set / = 2.5, p = 0.5/uJe,a = as before. 
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Fi gure 12: (a) An example of a trajectory that shows how the parametric resonance drive 
the reduced system from its almost equilibrium state to its open state, (b) A corresponding 
trajectory for the full model. 



For the initial condition, we set x a — qo(0) — 12.85, Uq — Po(Q) — 0,irj — 
0.86708 x u!q,(3q = and integrate the trajectory using the full equations, Eq. 



(29) where (<fc(0), P6(0)) is obtained by Eq. ( |47| ) and 9k(0) are determined by 
the Fourier modal transformation, Eq. ([3]). 

4.4.1 Remarks on This Class of Special Trajectories 

Here, we would like to make a few remarks: 



Despited its simplicity (when compared to the full system), the averaged 
reduced model is surprisingly accurate as illustrated by the fact that the 
three trajectories in Figure [TT^b) , Figure 12 'a), and Figure 12 'b) are very 



similar. Without a careful study of the averaged reduced equations, it 
may be difficult to guess that such a class of trajectoris will exist in the 
averaged reduced model, let alone in the reduced and the full model. 

It is also interesting to point out that the initial conditions for the reduced 
and the full models are essentially the same. Hence, if a small amount of 
initial energy is in a single Fourier mode, the dynamics of the reduced 
model of a two mode truncation looks very similar to the dynamics of the 



full equations of our DNA model. Figure 13 may illustrate this point in 
another way. This figure shows the projection of the same trajectory of 
the full model on the phase space of its first 15 modes. We observe that 

(i) the parametric excitation injects energy into the 6th mode, some of 
which transfers to the reactive mode and drive the full system from the 
almost equilibrium position to the region that marks the DNA division, 

(ii) only an extremely small amount of energy transfers from the excited 
mode to the other modes. This observation again shows that the two mode 
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Figure 13: Figure shows the projection of the same trajectory of the full model on the phase 
space of its first 15 modes. The parametric excitation injects energy into the 6th mode, some 
of which transfers to the reactive mode and drive the full system from the almost equilibrium 
position to the region which marks the DNA division. 

truncation can provide an adequate reduced model for studying the control 
of DNA division via parametric resonance. See Figure [4] for comparison. 

• Figure [l2| b) shows that the trajectory starts at the equilibrium position of 
the reactive mode but with certain small amount of energy in the trigger 
mode. Without the parametric excitation, the system will liberate near 
the equilibrium state if there is no friction or die down if the friction exists. 
However, if the parametric excitation is turned on at t = with the data 
provided above, the 1 : 2 parametric resonance will inject energy into the 
trigger mode, increase the value of /, make the trajectory to reach the 
region that marks the DNA division [x = 0) but with large energies in the 
trigger mode. 

Figure [14] shows a sequence of 5 snapshots of the evolution of this solution 
trajectory in the physical space when it is near the region that marks 
the DNA division. Notice that the DNA chain that corresponds to this 
solution is near its open state — - the chain is near the upright position 
(with its average angle near zero) but with a periodic swing. Figu res |14| 
(a) to (e) show one of its swings. Since the 6th mode is used, Figure [T4[a) 
shows a curve of pendula with six peaks (same for the other 4 figures). 



2N 



t=9344.7 



t=9345 








Figure 14: Figure shows a sequence of 5 snapshots of the evolution of the DNA chain near 
its opening state when its corresponding solution trajectory is near the region that marks the 
DNA division. 



Moreover, it is also interesting to point out that the time elapsed between 
Figure [T4|[ a) to Figure 14 e) is 2.7 units of time which is nothing but the 
period of the parametric excitation (27r/2a; with u> = 1.17). 

• Even though qualitatively the trajectories for the averaged reduced sys- 
tem and the reduced system look the same, quantitatively there are cer- 
tain discrepancy in their initial conditions. The main reason is that in 
computing the average of the Morse term M, the Taylor expansion at 
8 = is used. Hence, the equilibrium point for the reactive mode of 
the averaged reduced system is given by (x, y) = (12.59, 0) instead of 
(lOi Po) = (Vnd e ,0) = (12.85,0) for those of the reduced system. We ex- 
pect that if we compute M with another Taylor expansion at i/n9 e , the 
discrepancy will be much less. See §3.6[ However, since our concern at 



29 



this stage is to prove the concept that parametric resonance can be used 
to control the DNA division, we will not tackle this numerical issue for 
now. 

• For the cases where the initial energy is in more than one mode, say in the 
6th and the 7th modes, numerical simulation of the full model shows that 
this kind of trajectories still exist as long as one of the mode is dominant 
and the parametric excitation is in resonance with the dominant mode. 
This should not surprise us because while the parametric resonance will 
inject energy into the dominant mode, the friction will damp out the other 
mode. 

• More studies may be needed in the future for the tradeoffs between the 
amplitude /, the detuning parameter a, the frictional coefficient /i on one 
hand and the initial action-phase Iq,Po on the other. 

5 Conclusion 

In this paper we have studied the internal resonance, energy transfer, activation 
mechanism, and control of a model of DNA division via parametric resonance. 
Our study has been based on a methodology merging geometric reduction, par- 
tial averaging, techniques of chaotic transport, and control via parametric reso- 
nance. This methodology is not limited to our current model and its application 
can be extended to more general molecular and mechanical systems (possibly 
involving multiple scales). This study also highlights the importance of iner- 
tial effects in molecular dynamics, such effects are usually ignored in classical 
studies of molecular systems with Langevin equations (although Langevin equa- 
tions preserve the Gibbs distribution as an invariant distribution, they do not 
account for the electrostatic screening nor the hydrophobic effects of the solvent 
and they introduce strong and not necessarily justified assumptions on the dy- 
namic of molecular systems). This is why we have in this first study analyzed a 
noiseless system. Further studies are required to analyze the effects of inhomo- 
geneity, helicity (see [5]), and noise (we note again here that there is no unique 
way to introduce noise in such systems, henceforth the dynamical aspect of the 
model may become strongly biased without proper experimental validation). 
The inertial effects studied in this paper are important, not only from a general 
modeling aspect, but also because they can be targeted for purpose of control 
(possibly with low intensity electro-magnetic fields). There is also increasing 
evidence that these inertial effects play significant biological roles |39) . 

It would be interesting to extend our present framework to the models of 
DNA with helicity. As has been done in Ref . [51 [T7] , the effect of the helical 
geometry of DNA could be incorporated into the present model by introducing 
additional coupling between every N pendula (nucleotides), where N is typically 
4. This kind of coupling could induce another pathway for intramolecular energy 
transfer, which could in turn make the excitation of the reactive mode even more 
effective. The helicity of DNA could also be responsible for the coupling between 
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the dynamics of DNA bases and that of the DNA backbone. Since the division 
of DNA is associated with the slowest-scale dynamics among the dynamics of 
the bases as we have seen in the present study, the division dynamics of DNA 
could be coupled (or in resonance) with the slow-scale dynamics of the DNA 
backbone. In order to study this kind of couplings between the bases and the 
backbone of DNA, it would be important to implement a novel DNA model that 
takes into consideration the three-dimensional helical geometry more directly 
along the lines of Refs. [43l|36]. Since the helical geometry is ubiquitous among 
biomolecules, one can expect that the helicity plays a fundamental role in the 
functions of biomolecules. 
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A Partial Averaging of Lagrangian Equations is 
Equivalent to Partial Averaging of its Hamil- 
tonian for the Reduced Model 

A.l Partial Averaging of Lagrangian Equations 

Recall the reduced equations of motion in the Lagrangian form are given by 

q a +uj 2 a q a = -eM 1 {q ,q 1 ). (48) 
They can be rewritten as a first order system as follows: 
9o = V%>0 q 1 = Pj 

Po = -\feM Q Pl = -u*q 7 - eM 7 (49) 

Clearly, the set of equations in the reactive coordinates are already in the Stan- 
dard Lagrange Form with a small parameter y/e. Moreover, by using the angle- 
action variables defined by 



Pi — \J 2i 7 w 7 cos <pj (50) 

where </> 7 = uj-yt + ^ 7 , we transform the other set of equations also into the 
Lagrange Standard Form 

<jo = Vepo ^7= e(M 7 / y / 2i 7 w 7 ) sin <^ 7 

■po = —y/eMo 7 7 = — eM 7 J2/ 7 /w 7 cos</> 7 . (51) 



<7 7 = \ /2J 7 /w 7 sin 
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Hence, we can apply the standard theory of averaging (by averaging t from to 
27r/u; 7 ) and obtain the averaged reduced equaitons 



2tt 



r 2tr r2-K 



M # 7 /=— / cos <£ 7 d0 7 (52) 



y 2tt 7 ~""" 7 "' - 2tt 

where x = qo,y = po, I = I-f^ip — ^p 7 ,uj = lu^. 

A. 2 Partial Averaging of Its Hamiltonian 

Recall the averaged reduced Hamiltonian is given by 

ff 2 =~y 2 +u;I+~ jf M<% (53) 
where M(go,<?7) is a polynomial in qo,q 7 . Its Hamiltonian equation is given by 



" * =0J+ ilii u >'" 



It: 



e 

2ir dx 



d f 2 * 

Md<t>~, 1 = (54) 



A. 3 Two Methods are Equivalent for the Reduced Model 

Recall M = dM/dq , M 7 = dM/dq-,. Since 



dM _ dM dq 1 dM _ dM dq 1 

d(j) 7 dq 1 90 7 dl <9g 7 dl 



(55) 



and tp = 4> — u, the I and the ^ equations of Eq. (52) and Eq. (54 1 are the 
same. Moreover, since 



2?r a tix a /■2tt 



9M 9 
— — d0 7 = — . 
o dq dx J 



Md^y. (56) 



the two set of equations for the reactive coordinates are also the same if they 
are rewritten in the second order forms. 
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